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The quasiparticle spectrum of a two-dimensional d-wave superconductor in the mixed state, 
Hci <C -ff <C Hc2, is studied both analytically and numerically using the linearized Bogoliubov-de 
Gennes equation. We consider various values of the "anisotropy ratio" vf/va for the quasiparticle 
velocities at the Dirac points, and we examine the implications of symmetry. For a Bravais lattice of 
vortices, we find there is always an isolated energy-zero (Dirac point) at the center of the Brillouin 
zone, but for a non-Bravais lattice with two vortices per unit cell there is generally an energy gap. 
In both of these cases, the density of states should vanish at zero energy, in contrast with the 
semiclassical prediction of a constant density of states, though the latter may hold down to very low 
energies for large anisotropy ratios. This result is closely related to the particle-hole symmetry of 
the band structures in lattices with two vortices per unit cell. More complicated non-Bravais vortex 
lattice configurations with at least four vortices per unit cell can break the particle-hole symmetry 
of the linearized energy spectrum and lead to a finite density of states at zero energy. 



I. INTRODUCTION 



In the past few years several key experiments 
have demonstrated that the order parameter in high- 
temperature superconductors has d-wave symmetry 
rather than the conventional s— "vvave symmetry found in 
low-temperature superconductorsEJ. This fact has strong 
implications for the low-temperature thermodynamics of 
these systems, because the d-wave symmetry implies the 
existence of points on the Fermi surface where the gap 
function vanishes. At these points the bulk quasiparticle 
spectrum will be gapless and these states will be occupied 
even at very low temperature. 

High- Tc superconductors are extreme Type-II su- 
perconconductors and in a magnetic field larger than 
Hci develop a vortex lattice. The geometry of this 
vortex lattice haSr, been investigated via small ang|l|e 
neutron scatteringa, scajining tunneling microscopyS'cl 
and magnetic decorationO in YBa2Cu307 (YBCO) and 
BisSrsCaCuaOg (Bi2212). In YBCO twinned single crys- 
tals the vortex lattice in a magnetic field of 6 T parallel 
to the c axis looks like a skewed square lattice with an 
angle between primitive vectors of about 77°cl. Low-field 
magnetic decoration studies of Bi2212 between 70 G and 
120 G parallel to the c axisci find a vortex lattice very 
close to a hexagonal one. 

We have studied the quasiparticle spectrum of a d- 
wave superconductor in the vortex lattice within the 
framework of the Bogoliubov-de Gennes equation. The 
main question we have addressed is whether the spec- 
trum becomes gapped in the presence of a magnetic 
field Hci ^ ^ Hc2 and more generally what the 
energy spectrum looks like. This question was previ- 



ously approachfid via numerical simulatiouj-pf a tight- 
binding modeHu and semiclassical analysisO. Gor'kov 
and Schriefferd and, in a niore recent preprint using differ- 
ent arguments, Andersonll3 predicted that the quasipar- 
ticle spectrum of a d-wave superconductor in a magnetic 
field H <C Hc2 is characterized by broadened Landau 
levels with energy levels 

En=±hujHVn, n = 0,l,... (1) 

where luh ~ y2tJ^jAo/^. Here ujc — \eH\/mc is the 
cyclotron frequency and Aq is the maximum supercon- 
ducting gap. A key assumption of Gor'kov and Schrieffer 
and of Anderson, however, was to neglect the spatially 
dependent supcrfluid velocity which has-been shown to 
strongly mix Landau levels by Mel'nikovM. j— , 

Recently a preprint by Franz and Tesanoviola has given 
new insight into the problem. They introduced a gauge 
transformation that takes into account the supercurrent 
distribution and the magnetic field on an equal footing. 
In this way they map the original problem onto that of 
diagonalizing a Dirac Hamiltonian in an effective periodic 
vector and scalar potential with vanishing magnetic flux 
in the unit cell. Employing the Franz-Tesanovic trans- 
formation, we were able to tackle the problem of un- 
derstanding the band structure of this system, via both 
analytic and numerical methods. 

Our analysis will be limited to the spectrum at low en- 
ergies, in the case where the magnetic field is very small 
compared to 11^2 ■ In this case the distance between the 
vortices is large compared to their diameter, and we can 
ignore contributions to the spectrum from inside the vor- 
tex cores. The quasiparticle states of interest to us are 
then constructed from excitations close to the nodes of 
the energy gap of the zero-field spectrum, and we can 
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ignore mixing between different nodes. Tlius we ana- 
lyze tlie effects of tfie magnetic field and vortex lattice 
in a model where the zero field spectrum consists of four 
independent anisotropic Dirac cones. (We assume that 
the magnetic field itself is uniform in the sample, which 
is appropriate for a bulk superconductor provided that 
H ^ i/ci, but extends to even smaller fields for a very 
thin sample.) 

A very important numerical parameter in the prob- 
lem is the anisotropy ratio ao — vf/va- Here vp the 
Fermi velocity, while wa — ^o/pf is the quasiparticle ve- 
locity parallel to the Fermi surface, so the ratio ao mea- 
sures the anisotropy of the quasiparticle velocities at each 
Dirac point, in zero field. From angle-resolved photoe- 
mission spectroscopy and thermal conductivity measure- 
ments, the value of ao for higltTc superconductors turns 
out to be about 14 for YBCOy and 20 for Bi2212By. 
Conceptually, however, it is useful to consider the entire 
range of possible values for au , including the "isotropic 
case" = 1. This is particularly useful because we 
find that certain features of the spectrum, such as en- 
ergy gaps can become extremely small for large values of 
the anisotropy, and therefore become difficult to resolve 
numerically. 

Since each vortex in a superconductor carries only 
carry half of the normal flux quantum (f>o = ^c/|e|, it 
is necessary to choose a unit cell with an even number 
of vortices, so that the electron wavefunctions are single 
valued. Thus, if the vortices sit on a Bravais lattice with 
one vortex per unit cell, it is necessary to use a double 
unit cell, containing two vortices, in order to carry out 
the analysis. On the other hand, if the vortices sit on a 
non-Bravais lattice, with two vortices per unit cell, one 
can use directly the unit cell of the vortex lattice. 

Some key features of the quasiparticle band structure 
may be noted by looking at Figs. below, which give 
results of our numerical diagonalization (discussed in Sec- 
tion [v| ) for a square Bravais lattice, rotated by 45° from 
the quasiparticle anisotropy axis, with relatively small 
anisotropy ratio 1 < a/) < 4. One striking feature is 
the presence of band crossings both at zero and finite 
energy, regardless of the value of the anisotropy ratio 
ao- At least for relatively small anisotropy, the level 
crossings seem to be limited to the F and M point, as 
defined in Fig. |^. These band crossings are particularly 
interesting because in the absence of some special sym- 
metries, we would expect their probability to vanish for a 
two-dimensional Hamiltonian with broken time-reversal 



cally, at each point of the magnetic Brillouin zone, as 
can be seen from the plotted band structures and is fur- 
ther discussed in Section 0. Doi ng p erturbation theory 
calculations described in Section VII, we also found that 



symmetry, as is discussed in more detail in Section VIII 



We have studied the role played by various symme- 
tries of the Hamiltonian both to simplify our numerical 
computations and to try to understand why band cross- 
ings are allowed at some isolated points in the magnetic 
Brillouin zone. In particular, we focused on zero-energy 
states, as their existence changes qualitatively the ther- 
modynamic functions of the system, at very low temper- 
atures. For the lattices described above, we find exact 
particle-hole symmetry, both numerically and analyti- 



a crucial role is played by the Bravais nature of the vor- 
tex lattice. This symmetry, together with particle-hole 
symmetry is directly responsible for the spectrum staying 
gapless in the presence of a magnetic field. To prove this, 
we considered, both in perturbation theory and by nu- 
merical diagonalization, what happens if we deform the 
vortex lattice so that the distance between the A and B 
sublattices 2Rq defined in Fig. |l]is not 1/2 of the distance 
between two flux lines belonging to the same sublattice, 
measured along the diagonal. We found that particle- 
hole symmetry still exists at each point, separately, of 
the Brillouin zone but that, generally, gaps open up both 
at the F and M points, as can be seen in Fig. p. This 



result will be discussed in more detail in Section IX 



Another important question is whether there are any 
further zero-energy modes in addition to the ones at the F 
point. This is a very delicate question to address numeri- 
cally because it is hard to distinguish small gaps from real 
zero-energy eigenvalues, both because of finite numerical 
accuracy and, more importantly, because of finite grid- 
size effects. The latter can be particularly troublesome 
when dealing with lattice fermions as will be discussed at 
the beginning of Section |V^ . In numerical calculations, 
as noted by Franz and Tesanovic, it appears that for 
anisotropics of the order of 15, there are two lines in the 
Brillouin zone where the quasiparticle energy vanishes. 
Based on our symmetry analysis, however, we conclude 
that there actually remains a very small energy gap all 
along these lines, both for the case of a Bravais lattice 
and for a non-Bravais lattice with two vortices per unit 
cell. On the other hand, we find that for "rectangular" 
vortex lattices, isolated energy zeroes are allowed along 
the FX symmetry line (or -XF, by inversion symmetry). 

Different conclusions are reached if one allows for more 
complicated lattice structures, for example considering 
four vortices per unit cell. In this case it is possible to 
have superfluid velocity distributions without a center of 
inversion symmetry leading to non particle-hole symmet- 
ric energy spectra as shown in Fig. For large enough 
anisotropy, there is nothing that prevents lines of energy 
zeroes from appearing, and the density of states can be- 
come finite at zero energy. 

The scaling laws_for thermodynamic functions com- 
puted by Voloviklij witliin the semiclassical theory, 
and by Simon and Ledl3 in a framework closer to 
our appmach, have been tested experimentally (see for 
exampldlZl) . It is of interest to determine the relation be- 
tween the semiclassical approximation and the full quan- 
tum mechanical spectrum. We have extracted a den- 
sity of states from the band structures we computed for 
square vortex lattices and, contrary to the semiclassi- 
cal prediction of a constant density of states at zero en- 
ergy, we find a linearly vanishing density of states at low- 
energy for anisotropy ratio < 4, as shown in Figs. 
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Because of the previous discussion on the non-existence 
of Unes of energy zeroes, we are led to beheve that this 
result holds for any value of the anisotropy ratio as long 
as the vortex lattice has one or two vortices per unit cell. 
However the semiclassical approximation may be valid 
down to extremely low energies, when ao is large. 

In Section ^ we study the crossover between the semi- 
classical and quantum mechanical regions. There are two 
relevant energy scales which we call Ei and £2- For 
E > El the density of states is qualitatively identical to 
the bulk density of states in the absence of a magnetic 
field, although there are still noticeable features induced 
by van Hove singularities of the band structure in the vor- 
tex lattice. Below Ei, the presence of a magnetic field 
can be accounted for within a semiclassical approxima- 
tion, all the way down to an energy scale E2 where a full 
quantum mechanical calculation becomes necessary. 

We compare our numerically determined energy scales 
El and E2 witlL-the crossove«cales introduced by Kop- 
nin and VoloviUla'tj (see alsal3) Ei « fivp/d and E2^ w 
Tiv/^/d, where d is the average distance between vortices. 
While we agree with their expression for we notice 
that for large anisotropy our numerical analysis indi- 
cates that, at least for the geometries considered here, E2 
should go to zero much faster than l/ao- (The precise 
functional dependence of E2 on the anisotropy ratio ao 
will be the object of a paper currently in preparation.) 

The two energies Ei and i?^^ can also be writ- 
ten as El « Tcy/H/Hc2 and E^^ w Ei{Tc/Ef) = 
{T'^ / Ep'VJ H / Hc2- Experimentally, taking vp = 2.5 x 
lO^cm/f0 we find that Ei/iksVll) = 30 K/T^^^ and 
E^^/ikB^H) = 2K/T^/^ for both YBCO and,Bi2212, 
to a good approximation. Recent specific heatO and es- 
pecially low temperature thermal conductivitycJ experi- 
ments have been performed in the regime E < E^^ and 
still show good agreement with the semiclassical predic- 
tions, which also seems to suggest that the right crossover 
scale between the quantum mechanical and semiclassical 
regime is much smaller than E^^ , for large anisotropy 
ratios. Note that the Landau level energy scale (|l|) of 
Gor'kov and Schrieffer and of Anderson is approximately 
the geometric mean of Ei and E^^ . 

For E > El, the density of states is linear in energy 
with superimposed sharp peaks (logarithmic van Hove 
singularities) and the slope is the same as the one in 
zero magnetic field. When E2 < E < Ei, we are in 
the regime where the semiclassical theory predicts a con- 
stant density of states. This region shrinks as one low- 
ers the anisotropy ratio ao, and disappears entirely in 
the isotropic limit = 1 as is apparent by looking at 
Fig. ^. In the very low-energy limit E < E2, the semi- 
classical theory breaks down and the density of states has 
to be computed through the quantum mechanical spec- 
trum. We find that for E <^ E2 the density of states is 
linear and vanishes at zero energy, for the Bravais latitce. 

To summarize the structure of the paper, the model 
Hamiltonian is discussed in Sections |lTV. An analysis 



of particle-hole symmetry follows in Section In Sec- 
tion VI we describe the numerical methods and results of 



the band structure calculation for a Bravais vortex lat 
while in Sections 



tice, 



VII 



VIH we study the low-energy 
states in perturbation theory. More general vortex lat- 
tices with two and four vortices per unit cell are consid- 
ered in Section [X. In Section^, we focus on the density 
of states and compare it to the semiclassical predictions. 
Conclusions follow in Section Xl. 



II. LINEARIZED BOGOLIUBOV-DE GENNES 
EQUATION 

In a spatially inhomogeneous system, the standard ap- 
proach to the description of the quasiparticle spedxum is 
provided by the Bogoliubov-de Gennes equationEll. For 
an arbitrary gap operator (i.e. not necessarily s— wave), 
it reads TiedG^ = Eij} where ij) = {u, vf" is a Nambu 
2-spinor whose components arc the particlelike and hole- 
like part of the quasiparticle wave function, respectively. 
The B(>gpliubov-de Gennes operator TIbag we will con- 
sider ist3 



V 



2m 



-Ep 



A* 



A 

2m 



Ef 



(2) 



where Ep is the Fermi energy and m is the elec- 
tron effective mass. Notice that wa-are neglecting the 
self-consistent interaction potentialEJ (analogous to the 
Hartree-Fock potential in the normal phase) and any dis- 
order potential. The gap operator acts on components of 
the wave function as Ag{r) = J dr' Ad{r, r')g{r'). For a 
dxy superconductor, the gap operator can be expressed as 
A = ^{px,{py,A(r)}} where A(r) = Aoe'-^W and the 

brackets represent symmetrization {a, 6} — \{ab + ba). 
We choose this orientation instead of the more conven- 
tional dx2_y2 purely for notational simplicity; results do 
not depend on this choice. 

In the absence of a magnetic field, there are four points 
on the Fermi surface at p = {±pp,0) and p — {0,±pp) 
where the gap vanishes. If we are interested in the low- 
energy properties of the quasiparticle excitation spec- 
trum, we can linearize the Bogoliubov-de Gennes equa- 
tion around one of these points. This procedure is justi- 
fied because we are considering magnetic fields H <C Hc2 
and so the inverse magnetic length is much smaller than 
kp. We can choose to linearize around p = {0,pp) writ- 
ing the wave function ?/' — e^'^'^'^ip. The Bogoliubov-de 
Gennes equation will read {Him + 'Wrost)'0 — ^i^j with 
TYiin the leading linearized term 



_ vp{py-^^Ay) ^K,A(r)} 
«.in- 1 1 r„ -vp{py + ^^Ay) 



^fe,A*(r)} 



(3) 
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where vp = pf/tti is the Fermi velocity and Tircst is the 
remaining piece, which we wiU disregard being smaUer 

by O (yj^:^) ■ Notice that the hnearized eigenvalue prob- 
lem TiiinV' = Ell) is gauge covariant, so the linearized 
spectrum will be gauge invariant. 



III. FRANZ-TESANOVIC GAUGE 
TRANSFORMATION 

Following Franz and Tesanovi(S we can eliminate the 
phase factor e"^*^'^^ from the off-diagonal components of 
the Bogoliubov-de Gennes equation performing the sin- 
gular gauge transformation 
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FIG. 1. (a) A and B sublattices and vortex lattice unit cell, 
(b) The corresponding magnetic Brillouin zone. 
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gicl>Air) Q 







i4>B{r) 



(4) 



where (/)(r) = (/'yi(r) + (f>B{'r')- (/'yi(r) and (psir) are the 
contributions to the phase coming from the vortex sub- 
lattices A and B respectively, as defined in Fig. Q and 
will be computed later in the paper. 

The gauge transformed Bogoliubov-de Gennes opera- 
tor is 



VFPy VAPx 
VAPx -VFPy 
A 



VfV 



sy 



2 \^sx "sxJ 



2 (^sx Vgx) 



VFV, 



(5) 



where v\ = Aq/pf and the superfluid velocities corre- 
sponding to the A and B sublattices are defined as 



< = -(riV0,,--A), f,^A,B. 
m c 



(6) 



The operator (^) describes the dynamics of a free Dirac 
particle in a periodic potential. We can take advantage 
of the periodicity of the potential rewriting our spinors 
in Bloch form 



„ik r ( Uk 

Vk 



(7) 



where the functions Uk{r) and Vk{r) are themselves pe- 
riodic on the unit cell shown in Fig. |l] and the effective 
Hamiltonian acting on the Bloch spinors (U, V)"^ is 



n 



VFiPy 
VAiPx 



VJ 



hk 

A 



VfV 



sy 



va{Px + Tt-K) 

-VF{Py + TT-ky) 
2 \^sx "s: 



2 (^sx ^ss) 



VFVg 



(8) 



IV. ORDER PARAMETER SPATIAL 
DISTRIBUTION 

The last ingredient we need to determine the quasi- 
particle spectrum is the order parameter spatial distri- 
bution in the vortex lattice. The Bogoliubov-de Gennes 
method should in principle be used as a set of equations 
to be solved self-consistently, thus finding the quasiparti- 
cle spectrum and wavef unctions, and the spatially vary- 
ing order parameter distribution. As stated in the intro- 
duction, however, we restrict ourselves to the case where 
the vortex core size is small compared to the spacing be- 
tween vortices. Thus, we take the magnitude of the order 
parameter Aq to be constant everywhere except at the 
vortex sites, where it vanishes. We also assume that the 
magnetic field is constant in the sample. The phase (/)(r) 
of the gap function A(r) = Aoe*"^''"' is then obtained 
by minimizing a simplified Ginzburg-Landau free energy 
functionatd, which has the form 



const X 



(9) 



where we have taken e to be negative. The phase of the 



order parameter in the Landau gauge Ax 
is then given by 

zx {r -Ta) . 



-By, Ay 







= hm y . 



(10) 



where the sum runs over the vortex lattice. This series 
can be evaluated in closed form and the final result for 
the phase of the order parameter corresponding to the 
two sublattices (t)A{r) and (j>B{r) for a square vortex lat- 
tice with intervortex distance 2Rq ~ dy/2/2 is 



(/)A(r-) = Arg 



<t>B{r) 




(11) 



(12) 
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where QjJz,t) is the antisymmetric elhptic theta 
functioncS and the modular parameter t = i for a square 
lattice. These functions are not periodic, as they are 
gauge dependent, but using-thc properties of the theta 
functions under translationf£3 it is possible to show that 
the superfluid velocities vf'^{r) = m^^[h\/(j)A,Bi''') — 
{e/c)A] are. For a general lattice, with two vortices per 
unit cell, the Fourier representation of i'^(r) is given by 



Unk 
Vnk 



= S 



Vnk 



(15) 



<(r) = 



where Rq 



2TTh_ ^ iQ X g ,Q.(^„HM^ 
Q#o ^ 



fi^A,B (13) 



Rq and Rq — —Rq, Q are the reciprocal 



lattice vectors, and (P is the area of the magnetic unit 
cell. This in turn implies that the total superfluid ve- 
locity Vs = m~^[{h/2)V(j) — {e/c)A] is also a periodic 
function over the vortex lattice as is every gauge invari- 
ant quantity derived from it. 



V. PARTICLE-HOLE SYMMETRY 

Several symmetries play an important role in this prob- 
lem, both conceptually and computationally. Of course, 
the translational symmetry of the vortex lattice allowed 
us to introduce Bloch functions and recast the calcula- 
tion of the quasiparticle spectrum into a band theory 
framework. 

There are further symmetries that provide some help in 
understanding general features of the spectrum and sim- 
plify the diagonalization of the Bogoliubov-de Gennes 
equation (^). In the Introduction we mentioned that 
particle-hole symmetry and the Bravais nature of the vor- 
tex lattice are the two key ingredients that lead to the 
gaplessness of the quasiparticle spectrum. The impor- 
tance of th e Br avais lattice will be emphasized in detail 
in Sections VII -IX where the band structure is studied in 
perturbation theory. Here we want to focus on particle- 
hole symmetry. If (C/„fc(r), V„fc(r))^ is an eigenvector of 
the Bogoliubov-de Gennes equation (^) with eigenvalue 
Enk where n is a band index and A; is a wave vector in 
the first Brillouin zone, define 



Vnkir)^U:^{-r). 



(14) 



We claim that the spinor {Unk, Vnk)'^ is an eigenvector of 
the Bogoliubov-de Gennes operator (^) with eigenvalue 
—Enk- It is important to notice that in this way we are 
proving that particle-hole symmetry doesn't just hold on 
the whole spectrum as a set, it is an exact symmetry of 
the linearized Bogoliubov-de Gennes operator at every 
point in the Brillouin zone. In particular this means that 
the entire band structure should be exactly particle-hole 
symmetric. The proof of this statement is most easily 
constructed rewriting the transformation ([l^) more ex- 
plicitly in an operator notation 



with S = CTZT where C is the complex conjugation op- 
erator, TZ is the reflection through the origin operator 
and T = —ia2- It is easy to show that the linearized 
Bogoliubov-de Gennes operator defined in equation 
(^, and S anticommute {S^H} — using the property 
that the superfluid velocities for the two sublattices A 
and B are related by 



(16) 



in the coordinate system sketched in Fig. |l|, and vice 
versa. The particle-hole symmetry of the band structure 
is also observed explicitly in the numerical spectra as can 
be seen in Figs. for the case of a Bravais lattice. In 
the proof above, we only used the existence of a center 
of inversion in the vortex lattice. This holds also in the 
case of a non-Bravais lattice with two vortices per unit 
cell and so we expect to find a particle-hole symmetric 
band structure as well, which agrees with the numeri- 
cal results shown in Fig. ||. However, the proof fails, in 
general, for more complicated structures, with more than 
two vortices per unit cell where there is not a center of 
inversion symmetry anymore. 



VI. BAND STRUCTURE CALCULATIONS 
THROUGH NUMERICAL DIAGONALIZATION 

We have run extensive numerical diagonalization of the 
Hamiltonian (18), scanning the magnetic Brillouin zone 
for different values of the anisotropy ratio au = vf/va- 

Unlike Franz and Tesanovic, we decided to run our 
band structure calculations in real space, rather than mo- 
mentum space. The real space discretization leads to a 
sparse representation of the Hamiltonian which allows us 
to look at finer meshes then we would be able to in recip- 
rocal space. The algorithm we used for finding the eigen- 
functions and eigenvalues is a modified Lanczos-tyjpe 
method called the implicitly restarted Arnoldi methocO, 
implemented through the public domain Fortran 77 pack- 
age ARPACK. 

We note that the superfluid velocity Vg is singular at 
the vortex points (because we are taking the limit where 
the coherence length goes to zero), and its Fourier com- 
ponents (^2|) decay only algebraically. The real space 
and reciprocal space methods differ in the way in which 
they treat the large-wavevector cutoff. We think it is 
more direct to control the effects of this singularity in 
real space than in reciprocal space. However, we noticed 
a strong sensitivity of the results to the position of the 
vortex with respect to the real space mesh, until we cut 
the singularity off with a gaussian smoothing factor on a 
scale of a few grid points. Using a real space approach, 
the discretization of the superfluid velocity and it's reg- 
ularization near the vortices are controlled by separate 
parameters and can be optimized independently. 
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A disadvantage of using a real space approach for 
fermions is thaA one has to deal with the fermion dou- 
bling probleniEa. If one discretizes the Dirac Hamiltonian 
in two spatial dimensions in the most straightforward way 
using a rectangular mesh then one finds (in the absence 
of a scalar or vector potential) that there are 3 spuri- 
ous low-energy modes at the boundaries of the Brillouin 
zone of the mesh, i.e. when k/ir = (a~^, 0), (0, a^^) or 
(a~^, tty^), where and ay are the mesh spacings. This 
is a problem because the spurious modes get mixed with 
the physical low-energy modes (near fc = 0) by the inho- 
mogeneous potential. This problem has been thoroughly 
investigated in the lattice gauge theory community ajujl 
one way of getting rid of it was introduced by WilsonE3. 
Essentially, one introduces a fc— dependent mass term in 
the Dirac Hamiltonian that vanishes like fc^ at the center 
of the Brillouin zone and lifts the zero-energy modes at 
the boundaries of the unit cell to some high energy scale. 
Our choice for the Wilson term is 



, , 1 — cosfcrar , 1 — cosA:„a„ 
Hw = fiVF { K — + \ ^ ^ 



^2 



= hvF {XxO-xS^ + XyClySy) (72, 



(17) 



where 6^ and 6y are the second difference operators on the 
lattice. We are interested here in wavefunctions which 
vary smoothly on the scale of the vortex lattice spacing 
d. If we choose and a.y sufficiently small compared to 
d, for fixed Xx and Xy the Wilson term will have negligi- 
ble effect on the physical states, but the spurious states, 
which oscillate on the scale of a, will be pushed up to 
very high energies. 
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FIG. 2. Scaling analysis for the lowest three distinct eigen- 
values as a function of the Wilson parameter A in the isotropic 
case Of 15 ~ vf/va ~ 1 {Xx ~ ^y) for fixed number of mesh 
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When the superfluid velocities v'^{r) are included, the 
Wilson scheme breaks the symmetry that keeps the spec- 
trum gapless at the center of the Brillouin zone of the 
vortex lattice. Therefore, we have to perform a finite 
size scaling analysis to determine whether the small gaps 



we see in the numerics disappear in the limit A^:, Aj, 0. 
In Fig. 1^ one can see an example of such an analysis for 
the isotropic case ao = 1- If the spectrum is gapless 
one expects the Wilson term to open a gap linear in Xx 
and Aj,, which is very close to what we see for the lowest 
eigenvalue. We use this analysis to choose a value of A^; 
and Xy that is as small as possible but that doesn't get 
us in the region where we can see effects of the fermion 
doubling problem. 
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FIG. 3. Band structure and density of states for a square 
lattice with an = vf/va = 1. The dashed line is the spec- 
trum of the free Dirac Hamiltonian. Only the positive energy 
bands are plotted for clarity, the negative energy ones can be 
obtained by particle-hole symmetry (the density of states plot 
shows the overall particle-hole symmetry explicitly). Energies 
are in units of Ei = hvp/d. 

We have computed the band structure for several val- 
ues of the anisotropy ratio a]j. Figs. ||-^ show the band 
structure along symmetry lines and the density of states 
for the case of a square lattice with anisotropy = 1, 2 
and 4 respectively. The quasiparticle energy bands in a 
square vortex lattice are symmetric under the exchange 
of kx —kx or ky —ky and so only positive kx and ky 
have been considered. As we have already discussed in 
Section the bands are particle- hole symmetric, there- 
fore only positive energy bands are plotted. 

Noticeable features are the absence of a gap at the F 
point and further band crossings at higher energies also 
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at the r point. As we mentioned earlier and as will be 
analyzed in more detail in the following sections, the Bra- 
vais nature of the vortex lattice plays an essential role in 
keeping the spectrum gapless. The M point is also a spe- 
cial point, there are band crossings although not at zero 
energy. Also these band crossings will become avoided 
crossings once we modify the lattice into a non-Bravais 
one. 




FIG. 4. Band structure and density of states for a square 
lattice with oid = vf/va ~ 2. Negative energy bands can be 
obtained by particle-hole symmetry. 



VII. PERTURBATION THEORY AT THE F 
POINT 

One of the advantages of the Franz-Tesanovic gauge 
transformation is that it rephrases the problem in a form 
that is well suited to a perturbative analysis. The original 
linearized Bogoliubov-de Gennes equation doesn't easily 
separate into an exactly solvable unperturbed part plus a 
periodic (or quasi-periodic) perturbation while (^) is im- 
mediately recognizable as a two-dimensional free Dirac 
Hamiltonian perturbed by an effective periodic vector 
and scalar potential 

n = Ho + Hi (18) 
Ti-a = vpiPy + hky)a3 + ua(Px + 'hkx)(Ji (19) 

TYl V / A Fi \ (A R\ 

Wl = — \y,y -f V^y) (To + VF [Vsy - V^y) G 

+ ^'A {vt - ^^f.) ^i] (20) 



where Ci , i = 1 , . . . , 3 are the Pauli spin matrices and gq 
is the 2 by 2 identity matrix. 




FIG. 5. Band structure and density of states for a square 
lattice with a.u = vf/va —4. In figure (b), the dotted line is 
a smoothed out density of states on a scale AE — 0.25hvF /d, 
while the dashed line is the semiclassical density of states for 
a square lattice and od = 4. Negative energy bands can be 
obtained by particle-hole symmetry. 

The real space unit cell has to enclose an even number 
of vortices as each one of them carries only half a flux 
quantum $0 = hc/\e\. We will consider only unit cells 
with two vortices, but we will not restrict our analysis to 
Bravais lattices. The origin will be chosen at the center 
of inversion symmetry as shown in Fig. ^ and the po- 
sition of the two vortices Rq and —Rq is left arbitrary. 
The question we will try to address in perturbation the- 
ory is whether the presence of vortices introduces a gap 
in the quasiparticle spectrum. Then, if Tii is sufficiently 
small, we can limit our calculation to the center of the 
magnetic Brillouin zone (the T point). In Fourier space, 
the Hamiltonian ( |l^ ) at the F point can be written as 

+<-P.^3 + V«_p^ai (21) 
where the periodic potentials are 

vS'=2.^^cosg.ilo 
vS^=-2-^S-nQ.i?o (22) 



7 



V(3) -Ott^''^ 



sin Q ■ Rq 



The unperturbed eigenvalues are i?. 



v\h^Q'^ where Q 



(0) 

Q,± ~ 
nb2 are re- 



ciprocal lattice vectors. The basis vectors bi and 62 
are chosen such that, if ai and 02 are the basis vectors 
of the direct lattice, • bj = 2iTSij. The unperturbed 
eigenvectors can also be computed explicitly. If we write 
them in the form 



\E, 



(0) N 

■Q,±/ 



iQr 



(23) 



the coefficients a and (3 can be chosen real and have the 
symmetry property 



(24) 



and E^^^ 



and analogously for /3. 

We are interested in finding if the levels E_^_ 
are split by the perturbing potentials (^). Either even 
or odd orders in perturbation theory have to vanish be- 
cause if we change the sign of the superfluid velocity the 
splitting should remain unchanged as a consequence of 
particle-hole symmetry. It's easy to see that first order 
perturbation theory vanishes as the unperturbed eigen- 
vectors are constant spinors and the superfluid ve- 
locities v"^'^ have zero average. 



To find the effect of the perturbation to second or- 
der, we have to write an effective Hamiltonian for 
the two lowest energy bands. Using the formalism of 
Brillouin-Wigner perturbation theory, we may write the 
Schrodinger equation for a given wavevector in the form 



HcffiE)^ ^ E^ 



(25) 



where \1/ is a two-component spinor, and Ti.cS is a 2 by 2 
matrix defined by 

HcBiE) = Vk (ho + '^^ E-Ho^'Hi '^^ ^^^^ 

where Vk is the projection operator onto the two- 
dimensional subspace spanned by the unperturbed eigen- 
vectors Since we are looking at the behavior near 
zero energy, we can set = in (p6[). Also, to second 
order in Tii, we can neglect the Tii in the denomina- 
tor of the second term. Then, at the F point, where 
VqHqVo — (we are using the notation Vk=(o,o) ~ ^0 at 
the r point), we have 



n. 



(2) 
off 



-VoH: 



l-Vo 
Ho 



HiVo. 



(27) 



The matrix elements of this 2 by 2 matrix can be explic- 
itly calculated 



\-^o,+ l^ofr 1^0,+ 



(<il<l<l> 



E 



Q^O,i=± ^QA 



e: 



(0) 



(3) 



Q^Q,i = ± Eq ^ 



E 



(0) 



(3) 



V, 



(0) 

Q 



(28) 
(29) 
(30) 



Notice that for a Bravais lattice, Rq — (ai/4) 4- (02/4) 
and so if we write the reciprocal lattice vectors as Q = 
mbi + nb2, for m + n even Vq'' = Vq^ = while for 

m + n odd Vg'' = 0. Summing up the series in ( psf - 

|30| ) using the above mentioned property of the Fourier 

components of the potential and the symmetry of the 

unperturbed eigenfunctions ( p4[ ) it is easy to show that 
(2) 

Tipg vanishes, i.e. the perturbation does not open a gap 



in the spectrum up to second order. 

The third order contribution to the effective Hamilto- 
nian at the F point is 



(3) 



1 



K 



Vo 



Ho 



UiVo. 



(31) 



We will show that this matrix is also vanishing, to il- 
lustrate the point let's consider the off-diagonal matrix 
element 



E 



) + ^Q2-Qi("«i''i' 



2.«2 



J2A2 



(32) 
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Once again, we will use the Bravais lattice symme- 
try and consider pairs of terms in the sum correspond- 
ing to (Q,±) and (— Q,=f). In particular, if we write 
Qi 2 ~ '"^1,2^1 + and look at the relative sign of 

the terms in the sum (32) changing the signs of Qi and 
Q2 simultaneously, we find: 

Qi Q2 relative sign 

even even — 

even odd — 

odd even — 

odd odd — 

where "even" and "odd" refer to the parity of mi, 2 + 
ni_2- Adding up all the terms in pairs, we see that this 
matrix element vanishes as well as the diagonal ones, as 
can be easily checked. This calculation can be imme- 
diately extended to the fourth order, using exactly the 
same arguments and in fact, by induction, to any other 
order to show that to every order in perturbation theory 
the effective Hamiltonian at the center of the Brillouin 
zone for the lowest two bands vanishes. We thus find 
that to every order in perturbation theory the potential 
Til does not open a gap in the spectrum. 



VIII. PERTURBATION THEORY AWAY FROM 
THE r POINT 

The above analysis can be extended away from the 
r point. In particular we are interested in determining 
whether it's possible to find other points (possibly not 
symmetry points) in the Brillouin zone where the spec- 
trum is gapless. In the following, we will specialize our 
analysis to the case of a square lattice. pJBased on their 
numerical analysis, Franz and TesanovicE3 claim that for 
large enough anisotropy there is a whole line of zeroes 
that develops, in our notation, along a line parallel to the 
kx axis at a value of ky which depends on the anisotropy. 
(Note that our convention for the x and y axis is the op- 
posite of Franz and Tesanovic.) However, purely on sym- 
metry grounds, our effective Hamiltonian for the lowest 
two bands should be a complex herniitian 2 by 2 matrix. 
Particle-hole symmetry restricts the number of indepen- 
dent components to three (the effective Hamiltonian has 
to be traceless at every point in A:— space) but being in 
two dimensions we only have two parameters k^ and ky to 
vary. The system is obviously overdetermined and for a 
generic Hamiltonian of this kind we would not expect any 
zeroes, let alone lines of zeroes. The only way in which 
zeroes in the spectrum can develop is through some extra 
symmetry of the problem. We will see that there is such 
a symmetry only along the ky = axis. 

For a general wavevector k, at energy E = 0, we can 
write the effective Hamiltonian (Bq) in the form 



Heft - A{k)a3 + B{k)(7i + C(fc)cr2, 



(33) 



where cri, 172 and era are the Pauli spin matrices and A, B 
and C are real functions. 

In order for zero-energy states to exist A, B and C 
must all vanish simultaneously. We have seen that this 
happens at the F point, for a Bravais vortex lattice. To 
see if that can happen at other points, we will first con- 
sider the symmetry line kx = 0. We find that the coeffi- 
cient B{k) vanishes identically along this line. Although 
the coefficient A{k) is equal to vpky for the zeroth order 
Hamiltonian, it is possible that for large values of it 
could pass through zero and change sign at one or more 
values of ky other than ky = 0. If this occurs, and if C{k) 
were also zero, then there would be zero-energy states at 
these values of ky. However we shall see that along the 
line kx — 0, the coefRcient C{k) is different from zero 
in third order perturbation theory. Although C(fc) could 
have zeroes along the A:^: = axis for sufficiently large 
values of Qfu, there is no symmetry reason why these 
should occur at the points where A{k) vanishes. 

Similarly we find that along the ky — axis B{k) — 
C{k) = to all orders in perturbation theory, but that 
A(k) is generally non zero there. Isolated energy zeroes 
are therefore allowed by symmetry along the ky = axis 
and will be found if A{k) vanishes at any point on this 
symmetry line. 

For other points in the Brillouin zone, neither A nor 
B nor C vanish by symmetry, and there is no special re- 
lation between them. By varying kx and ky, one might 
find some isolated points where A and B vanish simul- 
taneously; however there is no reason why C should also 
vanish at such a point. Thus for a generic fixed value 
of ao, there should be no further zero-energy points in 
the Brillouin zone, other than along the ky = axis, 
where we find at least one state of zero energy at the F 
point. By varying an, however, it is possible that one 
could find special values where there are additional iso- 
lated zero-energy points. 

To summarize the results of the perturbative analy- 
sis, we find that there is always an energy zero at the 
F point, for a Bravais vortex lattice. In the case of a 
vortex lattice with a rectangular unit cell, rotated by 
45° from the quasiparticle anisotropy axis, there can be, 
for large enough anisotropy an, additional zero-energy 
states along the ky = axis. Of course, this result holds 
for any vortex lattice whose magnetic unit cell can be 
chosen as a rectangular unit cell properly oriented, in 
particular the triangular vortex lattice. At any other 
point of the magnetic Brillouin zone there will generally 
be no further zeroes in the energy spectrum, although 
there could be very low energy states. Also, at isolated 
values of the anisotropy ratio a d there could be energy 
zeroes at non-symmetry points in the Brillouin zone, but 
never lines of zero-energy states. 

We now show explicitly that for kx = the coefficient 
B{k) is zero to all orders in perturbation theory, while 
C{k) is nonzero at third order. We first consider the 
second order effective Hamiltonian 
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^ Ho / Pky,Q,i = /3(Q,,Q„+fej.j- Analogously to (|2§-^, we can 

where Vk is the projection operator onto the space calculate the matrix elements of this 2 by 2 matrix 



(4?.o.+iwrff^(^)i4:!o,+>=^, 



p(0) 
'(0,fc 

/(2)^Z7Mp(0) 



spanned by {|£^(o\),+>, l^fo, 



(0) 
(0,fc„ 



/(2) 



(0) 



E 



Let us define 



E 

1 



E-E^ 



(0) 



{EZ\jKsiEm\^)=E^ 



/(2) 



(0) 



X 

(o,fc„) - 



E-E 



(0) 



{aky,Q,zVj^^ +aky,Q,2Vj^^ +f3ky,Q,tV, 



{<^ky,QyQ^ - Pky.Q^Q - f3ky,Q,^V^Q^) 



E 



E-E 



(0) 



Pky,Q,i^Q' - f3ky,Qd^Q' + 0^ky,Q,i^Q 



(3) 



(1) 



(35) 



(36) 
(37) 



Looking at the off-diagonal terms ( p6| ) , we see that the 
imaginary part of these matrix elements vanishes if the 
vortex lattice is a Bravais lattice because, as we noted 



earlier, Vq^ vanishes when Vq' and Vq' do not and 
vice versa. Thus to second order, the coefficient of (T2 
is zero. This property holds for arbitrary fc, not just on 
the A;j,— axis, as afe,Q,i and Pk,Q,i are real numbers for 
every k and even though the off-diagonal matrix element 
away from the fc^, = axis has a more complicated struc- 
ture than in (Bq), its imaginary part will still be a sum of 



(1) 



;(3) 



polynomials in afc,Q,i and Pk.Q.i times Vq'^Vq'^\ which 
vanish identically. 

Going back to the kx — case, we can identify 



one further symmetry afe^,(„Q^,Qj,i = afe^,(Q,,Qj,j and 
l^ky,{~Q^.Qy).i = ~l3ky,{Q^.Qy).i wWch makcs the cti term 
in the effective Hamiltonian vanish -6(0, ky) = 0. 

If we look at the third order off-diagonal matrix ele- 
ment of the effective Hamiltonian, analogously to (|3^), 
we have, at i? = 0: 



\Eky,0,+ \Heff \Eky,0-) - 



(3)|p(0) 



E 

Qi,Q25^0,n,i2 = ± 



v(o) 



aky.Q 



2. '2 <'Q 



V, 



(1) 



E 



(0) 



E 



(0) 



X 



\^Ql-Qj^i'^ky,Qi,ii(^ky,Q2,i2 + Pky,Q^,iiPky,Q2A2) + '^Ql-Q^{^ky,Q^M^ky,Q2,i2 - Pky,Q^,hPky,Q2,i2) 

+ '^Ql-QS°^ky,Q^MPky,Q2a2 + Pky,Qi,h(^ky,Q2,i2] 



(38) 



Let us start by analyzing the coefficient of ai. Be- 
cause of the Bravais lattice symmetry, the matrix ele- 
ment ( |3^ ) is pure imaginary, and so can only contribute 
to CT2. This is true for every correction to the off-diagonal 
matrix element of the effective Hamiltonian coming from 
odd orders of perturbation theory: there will always be 
an odd power of V^*^^ Fourier coefficients in every term 
in the sum over intermediate states. The fourth order 
(and every other even order) correction could have a real 
part coming from terms with odd (meaning that if 
Q = {2t: / d){mx + ny), then to + n is odd) and every 
other reciprocal lattice vector alternating between even 
and odd. Considering pairs of terms of this kind with 
opposite Qx for every reciprocal lattice vector of the in- 
termediate states, it is easy to see that they will have 
opposite signs using the symmetry properties of the po- 
tential and of the unperturbed wave functions mentioned 
above. In particular, this implies that the ai term will 



keep vanishing along the fc^ = axis to every order in 
perturbation theory. 

Finally, we turn to the 12 term. This time, the Bra- 
vais lattice symmetry will ensure that corrections coming 
from even orders in perturbation theory are real numbers 
and so the only contributions to (J2 can come from odd 
orders in perturbation theory. The third order matrix 
element ( |3^ ) is the lowest non-vanishing one and it can 
be approximately evaluated numerically as a function of 
ky and the anisotropy ratio ao ■ We find that it is gener- 
ically non-zero when ky 0. 

The perturbative analysis for the ky — axis proceeds 
along the same lines. In this case, besides the Bravais 
nature of the vortex lattice, the symmetry responsible 
for the vanishing of the coefffcients B{k) and C(fc) is 

"fex,(Q^:-Qy):± = ~ Q^fe^c , (Q. , Q „ ) , T ^^'^ 1^ k, ,{Q ^ ,- Q y) ,± = 

/^fex,(Q;c.Qy).T' T^^'i matrix elements of the effective 
Hamiltonian can be explicitly evaluated taking Q!fc^_o,± — 
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±71 ^"^^ = 71-. 

To summarize, in this section we showed that for a 

Bravais lattice of vortices with a rectangular unit cell 
and for generic values of the anisotropy ratio zero- 
energy states can only be found along the ky = axis 
(one Dirac point is always found at the F point). Any- 
where else in the magnetic Brillouin zone we do not ex- 
pect to find further Dirac points, even though the gaps 
separating the particlelike and holelike bands could get 
very small. Also, for special values of it is not ruled 
out that there could be isolated energy zeroes anywhere 
in the Brillouin zone. 



IX. VORTEX LATTICE WITH A BASIS 

A. Two vortices per unit cell 

As we noted earlier, the key ingredient to find a gap- 
less spectrum at the center of the Brillouin zone in per- 
turbation theory was the Bravais nature of the vortex 
lattice. We can explore this connection further relaxing 
the constraint on the position of the vortices. The unit 
cell had to be doubled in order to enclose one quantum 
of magnetic flux $0 — hc/\e\; furthermore we divided 



the vortices into two sublattices A and B but kept them 
evenly spaced. We can leave the geometry of the two 
sublattices unchanged but displace them with respect to 
each other thus changing our lattice into a non-Bravais 
lattice with two vortices per unit cell. Let us assume that 
the two sublattices A and B are still square lattices with 
spacing d, but let us consider what happens if we let the 
distance between nearest A and B vortices be different 
from {^/2/2)d. For concreteness let us take 



■0 — 



1 - fx . 

: — X 



1-e. 



y d, 



as defined in Fig. 0. Then, Q Rq = (7r/2)[(l 
(1 — ey)'n\ and 



(3) 

-Q 



(39) 



c)m 



(40) 



The diagonal matrix elements of the effective Hamilto- 
nian at the F point (^8|) and (|3^), will keep vanishing 
because of the symmetries outlined above. On the other 
hand the off-diagonal terms do not vanish anymore and 
are purely imaginary, the second order effective Hamil- 
tonian is thus proportional to iT2- The series ( p9|) can 
be written as (keeping only the non-vanishing imaginary 
part) 



V, 



(0) 

Q 



E. 



(0) 



(1) 



E 

(m,n)5^(0,0),i=± 



(-1) 



ru+n+1 



2E, 



(0) 
Q,i 



■ sin7r(ea;TO -I- eyu) 



an 



(41) 



For small e we can expand the previous expression to 
first order in Crf. and Cy and using the symmetry 



^( — m,n) ,2 
/^( — m,n),2 ~ 



/^(m,n) ,i 



(42) 



it is easy to show that the series does not depend on ey. 
This result implies that, contrary to what we found in 
the Bravais lattice case, if the unit cell of the vortex lat- 
tice has a basis composed of two vortices with e^, ^ 0, 
the quasiparticle spectrum becomes gapped. In second 
order perturbation theory, the lowest eigenvalue at the 
center of the Brillouin zone depends linearly on ex for 
small distortions 



4% = ^9ic^D)\e.\ 
where the function g{aD) is defined as 



(43) 



giao) = vr- 



E 



(-1) 



(m>0,n),'i— ± 

1 



m 

Q,i 



'2oiQ.il3Q,im yfjQ j 

an 



(44) 



The asymmetry between and ey (present even in the 
isotropic = 1 case) is due to the term proportional 
to the identity matrix in the Hamiltonian (p^), in fact 

the terms in the series in ( pl| ) are proportional to VqK 
Linearizing the Bogoliubov~de Gennes equation around 
the Dirac points p = (pf,0) or p = (— ^^-,0) we would 
exchange the role of x and y so that for every distortion 
e the total density of states, defined as the sum of the 
density of states from the four Dirac points, exhibits the 
fourfold symmetry of the vortex lattice explicitly. 

In order to compare perturbation theory to the exact 
numerical diagonalization of the linearized Bogoliubov- 
de Gennes equation (nsh, let us consider the case e = 



ey. The matrix element ([4l|) can be evaluated nu- 



merically and for the case an 
and e = 0.1 we find 



1 (isotropic Dirac cone) 



while for 



n'^ = -0.055 



0.2 we have 



hvf 



-(72 



= -O.ll-f (72. 

a 



(45) 



(46) 
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The linear dependence of the gap on e is evident. 

We have run exact numerical diagonalization of the 
Hamiltonian (|l8|) and have found for the lowest energy 
eigenvalue in the e — 0.1 case 0.0765?!?;^/'^ while in 
the e = 0.2 case ^.X^lfivp jd which are in good agree- 
ment with the second order perturbation theory results. 
These gaps are larger than the gap induced by the Wil- 
son term in the band structure in a Bravais vortex lattice 
(where perturbation theory predicted gapless behavior) 
and scale linearly with e once the e = residual Wil- 
son gap is subtracted (for = 1 this is approximately 
Q.mShvpId). 

We also calculated the full band structure for the 
isotropic ao — 1, e = 0.2 case and the results are shown 
in Fig. 1^. Notice the gaps that open at the F and M 
points, while the rest of the band structure is qualita- 
tively unchanged, as one would expect from conventional 
perturbation theory. Contrary to the square lattice case, 
the band structure is not symmetric under the exchange 
of with —k^ or ky with 
energy bands in Fig. |[ 



ky as can be seen from the 




M Y r X M r -M 

FIG. 6. Band structure for a non-Bravais lattice with two 

vortices per unit cell in the ex = = 0.2, an = vf/va = 1 
case (solid curve). The dashed curve is the spectrum for a 
simple square lattice {ex = ey ~ 0). The gaps at the F and 
M point are real: they are much larger than the Wilson term 
contribution for anisotropy ao ~ 1, as can be seen compar- 
ing the dashed and solid curves. Notice that the symmetries 
kx ^ kx and ky * ky are broken in this case (the bands 
along the MF and -MF lines are different). Since particle-hole 
symmetry is still preserved at each point in the Brillouin zone, 
only the positive energy bands are plotted. 



B. Four vortices per unit cell 



In previous sections of the paper we have discussed 
the role played by particle-hole symmetry in determining 
some of the key features of the spectrum. In particular, 
we noticed that the vanishing of the density of states at 
zero energy, with or without a gap opening at the center 
of the Brillouin zone, is deeply related to this symme- 
try. To explore this connection further, it is of interest 
to find more complicated vortex lattice structures that 
break particle-hole symmetry and see the effect on the 
spectrum. 




FIG. 7. Band structure and density of states for a 
non-Bravais lattice with four vortices per unit cell with 
anisotropy ratio ao = vf/va = 5. The position of the vor- 
tices in the unit cell is displayed in the inset, notice that there 
is no center of inversion symmetry and thus particle hole sym- 
metry is absent. 



The density of states in the Bravais vortex lattice will 
be discussed in detail in the next section but it is clear 
that the gapless behavior at the F point implies the exis- 
tence of a low-energy window where the density of states 
is linear and vanishes at zero energy, just like in a homo- 
geneous d-wave superconductor. If the vortex lattice is 
distorted in the way discussed above, we find a gap in the 
quasiparticle excitation spectrum which depends on the 
magnitude and orientation of the distortion. The density 
of states will then vanish at zero energy in the general 
case of a vortex lattice with two vortices per unit cell. 



Particle-hole symmetry requires a center of inversion 
in the unit cell to exist. We can break this symmetry 
considering a unit cell with a basis consisting of four vor- 
tices, so that its area is 2d^. We choose a rectangular 
unit cell with sides —d < x < d, —d/2 < y < d/2. As an 
example, we can study the quasiparticle spectrum when 
the A-vortices are located at (0, 0) and {d/2, 0) and the B- 
vortices are located at (0, ztd/A), as shown in the inset in 
Fig. 0. In this case bands can go through the zero-energy 
axis away from cros sings or near-crossings and thus the 
discussion in Section VIII doesn't hold anymore and lines 



12 



of zeroes can be found. The band structure in Fig. ^ for 
anisotropy = 5 shows such hnes of zero-energy states. 
Instead of going to zero, the density of states stays fi- 
nite all the way down to zero energy, in close analogy to 
the prediction of the semiclassical theorj^. It may also 
be seen that in this example, there is no particle-hole 
symmetry for the overall density of states for the exhib- 
ited band structure. We must recall, however, that the 
exhibited states are derived from only one of the four 
Dirac points, p = (0,pi?), of the zero-field Fermi surface. 
If we include the contribution from the opposite point 
p = (0, ~pf) the overall particle-hole symmetry will be 
recovered. 




FIG. 8. Band structure and density of states for a 
non-Bravais lattice with four vortices per unit cell with 
anisotropy ratio q_d ~ vf/v/\ — 5. The position of the vor- 
tices in the unit cell is displayed in the inset, notice that 
there is a center of inversion symmetry with respect to which 
v^'^{—r) = —v^'^{r). Particle-hole symmetry is preserved, 
although not independently at every point in the Brillouin 
zone: states at fc and —k have to be exchanged. 

With four vortices per unit cell it is also possible to find 
superfluid velocity distributions that preserve particle- 
hole symmetry for the total density of states arising from 
a single Dirac point, but not at each point of the Brillouin 
zone separately as was the case for two vortices per unit 
cell. For example, the configuration depicted in Fig. ||, 
where the A-vortices are located at (±(i/2, 0) while the 
B-vortices are at the same position as before, shows such 
a distribution. Here, if we consider the transformation 



r —r, the superfiuid velocities v"^{r) and v^{r) are 
not exchanged as in equation (p^, rather they transform 
like v^'^{-r) = -v^'^{r). The Hamiltonian (||) goes 
into minus itself if we take r — > —r and simultaneously 
exchange k with —k. We then find E^{k) — —E^{—k), 
but no particle-hole symmetry at a fixed k. The density 
of states is an even function of the quasiparticle energy, 
as can be observed in Fig. ^. Lines of energy-zeroes are 
still allowed in this case, so the density of states may be 
finite at £^ = 0. 



X. DENSITY OF STATES: COMPARISON WITH 
THE SEMICLASSICAL THEORY 

The density of states is computed using a linear inter- 
polation for the band structure in between the sampled 
fc— points and is normalized as follows: 



:5{E-En{k)) 



(47) 



where the factor of 2 comes from spin degeneracy and n 
is a band index. As we noted in the previous section, this 
is the contribution to the total density of states coming 
from one of the four nodes of the zero-field quasiparticle 
spectrum. For a simple square vortex lattice in the orien- 
tation we are considering, the total density of states can 
be obtained simply by multiplying this result by four. In 
more complicated vortex lattices, a separate calculation 
of the density of states at one of the nodes rotated by 
90° with respect to p = (0,p_f) is generally necessary. 

If the vortex lattice is a Bravais lattice, the quasiparti- 
cle spectra are gapless regardless of the anisotropy ratio 
au and the density of states at very low-energy is linear 
in energy although the slope is renormalized by the po- 
tentials (and the renormahzation factor depends on «£>). 
The results of the numerical diagonalization are shown in 
Figs. I^ljfor ud = 1, 2, 4 respectively. The sharp peaks in 
the density of states are logarithmic van Hove singulari- 
ties: for topological reasons every band in two dimensions 
has at least two saddle points which contribute logarith- 
mically divergent peaks to the density of states. The van 
Hove singularities show up as finite-height peaks in the 
numerical evaluation of the density of states because of 
the linear interpolation scheme used for the band struc- 
ture. Averaging these peaks out, however, one can see 
that at high energy the density of states reproduces the 
behavior expected for the quasiparticles in the absence 
of a magnetic field N{E) — \E\/ (nh'^vpVA) as shown in 
Fig. H for the = 4 case. 

We want to compare our results fj?-fhe semiclassical 
picture studied primarily by Volovikot3. This approach 
takes into account the superfluid velocity Vg distribution 
through the Doppler shift of the quasiparticle energy 



E{k, r) = ±Jek + Ad(fc)2 + hk ■ v,{r), 



(48) 
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where is the kinetic energy measured with respect 
to the Fermi surface. Within this framework KopHiB, 
and Volovik introduced two crossover energy scaleal^jla'ta 
El « hvp/d and « TiVA/d. The first energy scale 

El « hvp/d marks the boundary between the temper- 
ature dominated regime and the superflow dominated 
regime in the thermodynamic functions. Physically, this 
crossover corresponds to the WKB eigenfunctions becom- 
ing extended on a scale comparable to the intervortex dis- 
tance, at least in one direction. For E > Ei, the states 
are unaffected by the magnetic field and the density of 
states is linear with the same slope one would find in 
the bulk without a vortex lattice, as we discussed in the 
previous paragraph. For E < Ei, the semiclassical den- 
sity of states is essentially independent of energyO and 
for our order parameter distribution on a square lattice 
we calculate it to be 

Volovik finds essentially the same resulti, with an unde- 
termined numerical prefactor which depends |-©n the ge- 
ometry of the vortex lattice. Won and MakiE3, using a 
somewhat different model, calculate this geometric pref- 
actor for a general vortex lattice structure (although only 
Bravais lattices are considered) and find a result of the 
same order of magnitude of (Eoh, for a square lattice. 

The second crossover marks the boundary where a 
full quantum mechanical picture becomes important and 
the semiclassical analysis breaks down. We will call 
this scak—Ea- As we mentioned above, Kopnin and 
Volovikta'tSllj argue that this scale is linear in l/ao, 
and is given by ~ hv^/d and, for energies in the 

range E2^ < E < Ei, they predict a constant density 
of states. In terms of band structure this means that we 
need to find a direction in fc-space in which the bands 
are flat on a scale of i?2- If we assume that the pertur- 
bation induced by the magnetic field is weak, and thus 
that the band structure is only weakly renormalized, we 
find that for large enough anisotropy 3> 1 several 
bands will start overlapping for energies E < Ei and will 
have a small dispersion in the kx direction of the order 
of hvA/d. The density of states for E < Ei would then 
be essentially constant down to energies of the order of 
hvA/d, which corresponds exactly to E^^ . For energies 
E < E2^ we would have just a single band and the 
density of states would drop linearly to zero. The key as- 
sumptions that enter in this argument (weak perturbing 
potential) are essentially the absence of vanishing ener- 
gies, or nearly vanishing energies, at any point in the 
Brillouin zone other than the F point and a small renor- 
malization of the slope of the lowest energy bands at the 
F point in any direction of the Brillouin zone. These as- 
sumptions are satisfied for rather small anisotropy ratios 
an <C 10 but seem to fail for higher values of a^i, as we 
will show. 

In Fig. we can see how the semiclassical description 
starts developing. For a_D = 1 there is no resemblance to 



the semiclassical behavior for any energy, the bands are 
very distinct at low energy and there is only one crossover 
from a quantum mechanical region to a purely classical 
region (E > Ei) without any hint of constant density of 
states. Changing the anisotropy to a^i = 2 or even better 
ao = 4 a hint of the semiclassical region starts opening 
up and one can identify a trend towards a flat density of 
states between roug hly E^^ and Ei. For energies much 
lower than i?!"^ the density of states is linear and goes 
to zero at zero energy. Although there is only one band 
(even in the ao = 4 case) below Ei , the semiclassical de- 
scription seems to start working remarkably well. From 
these plots we can already notice discrepancies with the 
Kopnin and Volovik picture. Notice that, while for the 
ao = 2 case the ratio between the slopes in the ky and 
kx directions is Vp / 2.5 ^ an (the superscript R in- 
dicates that these are the renormalized velocities in the 
two above mentioned directions in fc-space), already in 
the ao — 4 case the same ratio is roughly 17 which is 
much larger than ud- The scale at which the flat den- 
sity of states should break down, E2, seems to be much 
smaller than the simple argument above would predict. 
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FIG. 9. Band structure for a square lattice with anisotropy 
ratio an = vf/va = 8, 12, 15 (solid, dashed and dot-dashed 
lines respectively) along the symmetry line FY. Notice the ap- 
parent vanishing of the energy spectrum at a wavevector close 
to the Y point, for the large values of an. (Our numerical 
calculation cannot distinguish between a true zero and a very 
small energy gap.) A second small energy gap is developing 
about one third of the way from F to Y. Only the positive 
energy bands are plotted for clarity, the negative energy ones 
can be obtained by particle-hole symmetry. In the an ~ 15 
only the two lowest bands are plotted, while all the bands 
with energy E < Tivp/d are plotted for an = 8, 12. The gaps 
at the r point are fictitious and are due to the Wilson term 
in the Hamiltonian. 

Besides a large renormalization of the slope of the en- 
ergy bands at the F point, which occurs even for rel- 
atively small anisotropics, we find for large anisotropics 
that there are lines in the Brillouin zone where the energy 
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of the lowest ban d is v ery close to vanishing. As we dis- 
cussed in Section VII] , we do not expect to find points, 
other than the F point and possibly along the ky = 
axis, where the energy is exactly zero (except, conceiv- 
ably, at isolated values of the anisotropy ratio ao)- For 
anisotropies > 8 we cannot resolve numerically any 
dispersion along the direction, so only the FY line of 
the band structure carries information. In Fig. ^ we have 
plotted a few of the lowest energy bands corresponding to 
values of the anisotropy ratio ao =8, 12, 15. One can im- 
mediately see one of these energy near-zeroes developing 
along the FY axis. In general, the crossover scale E2 will 
be set by the larger of the energy gap at this point and 
the energy dispersion in the direction. The density of 
states will be very close to a constant for energies larger 
than E2 and will drop towards zero with decreasing en- 
ergy for _B <C i?2- In this way the high- anisotropy limit 
approaches the semiclassical prediction much faster then 
linearly in l/ao- (The precise functional dependence of 
E2 on an will be discussed in a later publication.) 

For more complicated lattices, where there is not 
particle-hole symmetry at each point in the Brillouin 
zone, there is no argument to prevent zero crossings, and 
we do indeed find lines of zeroes for large values oi Ud 
(see Figs. |^ and ^. Thus the density of states is finite 
at zero energy and the semiclassical results may apply 
down to zero energy. 



XI. CONCLUSIONS 

In conclusion, we have studied the quasiparticle spec- 
trum of a d-wave superconductor in the mixed state. 
One important step in solving the problem has-, been 
the transformation due to Franz and Tesanovicll3 that 
maps the original Bogoliubov-de Gennes equation into 
a Dirac Hamiltonian in an effective periodic vector and 
scalar potential corresponding to zero average magnetic 
field. We have found both numerically and in perturba- 
tion theory that for a Bravais lattice of vortices the spec- 
trum remains gapless when a magnetic field is turned on. 
We have showed that for vortex lattices which preserve 
the particle-hole symmetry of the energy spectrum there 
can only be other isolated Dirac points in the Brillouin 
zone and so they cannot change qualitatively the very 
low-energy density of states. Different conclusions are 
reached when more complicated vortex lattice structure 
(for example with four vortices per unit cell) are consid- 
ered where lines of energy-zeroes can be found. In this 
case, the density of states is finite at zero energy for large 
enough anisotropy ratio . A non-Bravais vortex lattice 
with two-vortices per unit cell can break the symmetry 
that keeps the spectrum gapless and open gaps whose 
magnitude depends on both the magnitude and orien- 
tation of the distortion observable both in the numerics 
and perturbation theory, with good agreement between 
the two. Finally, the high-anisotropy limit has been in- 



vestigated and it's relation to the semiclassical analysis 
explained. The crossover scale between the semiclassi- 
cal and quantum mechanical regime E2 goes to zero for 
large values of the anisotropy much faster than linearly 
in l/ao, at least for the vortex lattice geometries consid- 
ered here, and the density of states quickly approaches a 
constant value for energies E2 < E < Ei. 
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